%% Directory

% Matlab needs to be configured such that it can find the fun_***.m files
% one way to do this is 'addpath'
addpath 'LOCATION OF THE MATLAB .m FILE FOLDER HERE'


%% Exogenous parameters
clear all

T=85;
t=58;
t_retstart = 62; % This is when income starts to decline near retirement
t_ret = 66;
d=2; % dimensions
beta=0.98;
R_base = 1/beta;
R_shocked = R_base - 0.01;

logshifter = 10/1000000; % as in empirical spec

w = NaN(T,1);
p = NaN(T,1);

s_vector = (1:1:T)';

p_base = NaN(T,1);

grossinc = 0.77;


netinc = grossinc*(1-0.40); % net of assumed income tax rate
initialwealth = 4.23;

p_base(t:t+5) = [0.72;0.73;0.73;0.75;0.76;0.77];
p_base(t+6:T) = 0.77;

w(t) = netinc + initialwealth;
w(t+1:t_retstart) = netinc;
aidur = 0.5*1/(t_ret-t_retstart); % annual income decrease until retirement
w(t_retstart:T) = netinc - min(netinc*aidur*( s_vector(t_retstart:T)-t_retstart), netinc*0.5);

c_hat = 0.93*grossinc*(1-0.40);
g_hat = 5225*0.30/1000000;


%% Basline (No-Entry-Cost) Model. Giving Response for Different Values of the EIS
varepsilon_vector = [0.09 0.44 0.88]'; % Approximately 95% CIs around point estimate
sigma_vector = [(0.01:0.005:0.06)';0.061;0.062;0.063;0.06303;0.06305;0.064;0.065;0.066;0.067;0.068;0.069;0.07;0.0701;0.0702;0.0703;0.0704; 0.0705;0.0707;0.071;0.072;0.073;0.074;0.075;0.077;0.080;(0.0810:0.0001:0.082)';(0.0830:0.001:0.089)';0.09;(0.1:0.01:0.40)';(0.42:0.02:0.8)';(0.8:0.1:1.1)'];

p = p_base;
Dlog_g = NaN(size(sigma_vector,1),size(varepsilon_vector,1));

for ivarepsilon = 1:size(varepsilon_vector,1)
    varepsilon = varepsilon_vector(ivarepsilon);
    for isigma=1:size(sigma_vector,1)
        sigma = sigma_vector(isigma);

        kappa = c_hat^(-1/sigma)*g_hat^(1/varepsilon)*p(t);
        
        % Solve for unshocked g
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        R = R_shocked;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

        % Calculate treatment effect
        Dlog_g(isigma,ivarepsilon) = log(logshifter+gt_shocked)-log(logshifter+gt_base) 

        if isinf(kappa)
            Dlog_g(isigma,ivarepsilon) = NaN;
        end
    end
end

%% Export table to stata to make figure

[[NaN; sigma_vector],[varepsilon_vector'; Dlog_g]]

% This is manually pasted into wdon1_calibration_stata_plot.do, 
% below the line that only consists of 
% '	0		.			.			.		.			.'
% and before the next line that says 'end'


%% Calculate implied UNCOMPENSATED own-price elasticity
p_shocked = NaN(T,1);
varepsilon = 0.44;
sigma =0.0815;

% shock is 10% price increase
p_old = 0.73;
p_new = 1;
p_base(t:T) = p_old;
p_shocked(t:T) = p_new;


% BASELINE uncompensated elasticity
kappa = c_hat^(-1/sigma)*(g_hat)^(1/varepsilon)*p(t);
        % Solve for unshocked g
        p = p_base;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        p = p_shocked;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

        % Calculate uncompensated own-price elasticity
        e=(log(gt_shocked)-log(gt_base))/(log(p_new) - log(p_old))


% INCREASING BASELINE AMOUNT OF GIVING --> uncompensated elasticity
kappa = c_hat^(-1/sigma)*(g_hat*10)^(1/varepsilon)*p(t);
        % Solve for unshocked g
        p = p_base;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        p = p_shocked;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

        
 % Calculate uncompensated own-price elasticity
 e=(log(gt_shocked)-log(gt_base))/(log(p_new) - log(p_old))


%% Calculate Effect of Wealth Tax on Uncompensated Elasticity
p_shocked = NaN(T,1);
varepsilon = 0.44;
sigma =0.0815;

% own-price shock is 10% price increase
p_old = 0.73;
p_new = 1;
p_base(t:T) = p_old;
p_shocked(t:T) = p_new;

% R shock is 1% linear wealth tax
R_shocked = R_base - 0.01;

% BASELINE uncompensated elasticity (no wtax)
kappa = c_hat^(-1/sigma)*(g_hat)^(1/varepsilon)*p(t);
        % Solve for unshocked g
        p = p_base;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        p = p_shocked;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

% Calculate uncompensated own-price elasticity
e=(log(gt_shocked)-log(gt_base))/(log(p_new) - log(p_old))




%  Uncompensated elasticity with proportional wealth tax
        % Solve for unshocked g
        p = p_base;
        R = R_shocked;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        p = p_shocked;
        R = R_shocked;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

% Calculate uncompensated own-price elasticity
e=(log(gt_shocked)-log(gt_base))/(log(p_new) - log(p_old))


%% MPG out of wealth taxes 


p_base(t:t+5) = [0.72;0.73;0.73;0.75;0.76;0.77];
p_base(t+6:T) = 0.77;
sigma = 0.086;
varepsilon=0.44;

kappa = c_hat^(-1/sigma)*(g_hat*1)^(1/varepsilon)*p(t);

  % Solve for unshocked g
        p = p_base;
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        p = p_base;
        R = R_shocked;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

        % Also need to calculate denominator: effect on wtax, use FOCs and
        % budget constraint
        ct_shocked = (p(t)/kappa)^sigma * gt_shocked^(sigma/varepsilon);
        stplus1 = w(t)-p(t)*gt_shocked-ct_shocked;
        wtax = -stplus1*(R_shocked-R_base);
 
 % MPG out of wealth taxes       
 MPG = (gt_shocked-gt_base)/wtax

%% Create data for plot of how MPG varies with EIS for different Kappas


p_base(t:t+5) = [0.72;0.73;0.73;0.75;0.76;0.77];
p_base(t+6:T) = 0.77;

varepsilon = 0.44;
sigma_vector = [(0.02:0.0025:0.0699)';0.07;0.0705;0.071;0.075;(0.076:0.001:0.081)';(0.0811:0.0001:0.08149)';(0.085:0.001:0.25)';(0.26:0.01:0.40)';(0.42:0.02:0.8)';(0.8:0.1:1.5)'];


p = p_base;

Nkappa = 2;


MPG_vector= NaN(size(sigma_vector,1),Nkappa);

for ikappa = 1:Nkappa
    for isigma=1:size(sigma_vector,1)
        sigma = sigma_vector(isigma);
        
        % set kappa
        if      ikappa == 1
            kappa = c_hat^(-1/sigma)*g_hat^(1/varepsilon)*p(t); % original kappa
        elseif  ikappa == 2
            kappa = (c_hat-0.05*grossinc)^(-1/sigma)*(0.05*grossinc)^(1/varepsilon); % 5% giving kappa
        end

        % Solve for unshocked g
        R = R_base;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);

         % Solve for shocked g
        R = R_shocked;
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );

        syms x
        gt_shocked = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);
        
        % calculate wtax payments of treated 
        ct_shocked = (p(t)/kappa)^sigma * gt_shocked^(sigma/varepsilon) ;
        stplus1 = w(t)-p(t)*gt_shocked-ct_shocked;
        wtax = -stplus1*(R_shocked-R_base);
        
        MPG_vector(isigma,ikappa) = (gt_shocked-gt_base)/wtax;

    end
end

[sigma_vector MPG_vector]

% This is manually pasted into wdon1_calibration_stata_plot.do, 
% below the line that only consists of '0.09 . .'
% and before the next line that says 'end'

    